#include "SurfaceExtractor.h"
#include <fstream>
#include <sstream>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <iomanip>

SurfaceExtractor::SurfaceExtractor(const GridModel& model)
    : gridModel(model) {
}

bool SurfaceExtractor::isZCoordEqual(double z1, double z2) const {
    return std::abs(z1 - z2) < zTolerance;
}

bool SurfaceExtractor::isEdgeCoincident(const std::vector<double>& edge1Zs, 
                                      const std::vector<double>& edge2Zs) const {
    if (edge1Zs.size() != edge2Zs.size()) return false;
    for (size_t i = 0; i < edge1Zs.size(); ++i) {
        if (!isZCoordEqual(edge1Zs[i], edge2Zs[i])) {
            return false;
        }
    }
    return true;
}

std::vector<double> SurfaceExtractor::getPillarZCoords(int i, int j) const {
    const auto& coordData = gridModel.getCoordData();
    const auto& gridSize = gridModel.getGridSize();
    
    int idx = (i + j * (gridSize.Nx + 1)) * 6;
    return {coordData[idx + 2], coordData[idx + 5]};
}

void SurfaceExtractor::loadFaultData(const std::string& filename) {
    std::ofstream logFile("fault_data_log.txt");
    if (!logFile) {
        throw std::runtime_error("无法创建断层数据日志文件");
    }

    std::ifstream file(filename);
    if (!file) {
        logFile << "错误: 无法打开断层文件: " << filename << std::endl;
        logFile.close();
        throw std::runtime_error("无法打开断层文件: " + filename);
    }

    logFile << "开始读取断层数据文件: " << filename << std::endl;
    logFile << "文件是否打开: " << (file.is_open() ? "是" : "否") << std::endl;
    logFile << "文件是否正常: " << (file.good() ? "是" : "否") << std::endl;

    // 记录MAPAXES参数
    const auto& mapAxes = gridModel.getMapAxes();
    logFile << "\nMAPAXES参数：" << std::endl;
    logFile << std::scientific << std::setprecision(6);
    logFile << "Y轴点: (" << mapAxes.yAxisX << ", " << mapAxes.yAxisY << ")" << std::endl;
    logFile << "原点: (" << mapAxes.originX << ", " << mapAxes.originY << ")" << std::endl;
    logFile << "X轴点: (" << mapAxes.xAxisX << ", " << mapAxes.xAxisY << ")" << std::endl;

    std::string line;
    FaultData currentFault;
    bool readingFault = false;

    while (std::getline(file, line)) {
        if (line.empty()) continue;

        std::istringstream iss(line);
        std::string firstWord;
        iss >> firstWord;

        if (firstWord == "FAULT") {
            if (readingFault && !currentFault.nodes.empty()) {
                faultData.push_back(currentFault);
                logFile << "保存断层 " << currentFault.name << "，共 " 
                       << currentFault.nodes.size() << " 个节点" << std::endl;
            }
            currentFault = FaultData();
            std::string name;
            iss >> name;
            currentFault.name = name.substr(1, name.length() - 2);
            readingFault = true;
            logFile << "\n开始读取新断层: " << currentFault.name << std::endl;
        } 
        else if (readingFault) {
            std::string line_copy = line;
            std::istringstream data_iss(line_copy);
            double x, y, column_id, z;
            
            if (data_iss >> x >> y >> column_id >> z) {
                logFile << "\n处理节点数据：" << line << std::endl;
                logFile << "原始大地坐标: (" << std::fixed << std::setprecision(6) 
                       << x << ", " << y << ", " << z 
                       << "), 柱线ID: " << static_cast<int>(column_id) << std::endl;

                // 计算相对于原点的偏移
                double deltaX = x - mapAxes.originX;
                double deltaY = y - mapAxes.originY;
                logFile << "相对原点偏移: (" << deltaX << ", " << deltaY << ")" << std::endl;

                // 计算变换矩阵
                double dx = mapAxes.xAxisX - mapAxes.originX;
                double dy = mapAxes.xAxisY - mapAxes.originY;
                double length = std::sqrt(dx * dx + dy * dy);
                double ex = dx / length;
                double ey = dy / length;

                // 转换坐标
                auto transformed = gridModel.transformCoordinates(x, y);
                
                FaultNode node;
                node.x = transformed.first;
                node.y = -transformed.second;  // 对y坐标取负值
                node.z = -z;  // z坐标仍然取负值
                node.column_id = static_cast<int>(column_id);
                node.fault_name = currentFault.name;

                logFile << "转换后坐标: (" << node.x << ", " << node.y << ", " << node.z 
                       << "), 柱线ID: " << node.column_id << std::endl;

                currentFault.nodes.push_back(node);
            }
        }
    }

    // 保存最后一个断层
    if (readingFault && !currentFault.nodes.empty()) {
        faultData.push_back(currentFault);
        logFile << "\n保存最后一个断层 " << currentFault.name << "，共 " 
               << currentFault.nodes.size() << " 个节点" << std::endl;
    }

    logFile << "\n断层数据统计：" << std::endl;
    logFile << "总断层数: " << faultData.size() << std::endl;
    for (const auto& fault : faultData) {
        logFile << "断层 " << fault.name << ": " << fault.nodes.size() << " 个节点" << std::endl;
        for (const auto& node : fault.nodes) {
            logFile << "  柱线 " << node.column_id << ": (" 
                   << node.x << ", " << node.y << ", " << node.z << ")" << std::endl;
        }
    }

    file.close();
    logFile.close();
}

void SurfaceExtractor::createFaultSurfaces() {
    std::ofstream logFile("fault_surfaces_log.txt");
    if (!logFile) {
        throw std::runtime_error("无法创建断层面片日志文件");
    }

    logFile << "开始创建断层面片" << std::endl;
    logFile << "断层数量: " << faultData.size() << std::endl;

    faultQuads.clear();

    for (size_t faultIndex = 0; faultIndex < faultData.size(); ++faultIndex) {
        const auto& fault = faultData[faultIndex];
        logFile << "\n处理断层 " << fault.name << " (索引: " << faultIndex + 1 << ")" << std::endl;
        
        // 按柱线编号分组并排序节点
        std::map<int, std::vector<FaultNode>> columnNodes;
        for (const auto& node : fault.nodes) {
            columnNodes[node.column_id].push_back(node);
        }
        
        logFile << "柱线数量: " << columnNodes.size() << std::endl;
        
        // 对每个柱线内的节点按Z坐标排序
        for (auto& pair : columnNodes) {
            std::sort(pair.second.begin(), pair.second.end(),
                [](const FaultNode& a, const FaultNode& b) {
                    return a.z < b.z;
                });
            logFile << "柱线 " << pair.first << " 包含 " << pair.second.size() 
                   << " 个节点" << std::endl;
        }

        // 创建断层面片
        std::vector<int> columnIds;
        for (const auto& pair : columnNodes) {
            columnIds.push_back(pair.first);
        }
        std::sort(columnIds.begin(), columnIds.end());

        int quadCount = 0;
        // 在相邻柱线之间创建四边形
        for (size_t i = 0; i < columnIds.size() - 1; ++i) {
            const auto& currentColumn = columnNodes[columnIds[i]];
            const auto& nextColumn = columnNodes[columnIds[i + 1]];

            size_t numNodes = std::min(currentColumn.size(), nextColumn.size());
            
            logFile << "处理柱线对 " << columnIds[i] << " - " << columnIds[i + 1] 
                   << ", 可创建 " << (numNodes - 1) << " 个四边形" << std::endl;
            
            for (size_t j = 0; j < numNodes - 1; ++j) {
                if (faultIndex >= 9) {
                    logFile << "跳过断层索引 >= 9 的面片" << std::endl;
                    continue;
                }
                
                Point3D p1(currentColumn[j].x, currentColumn[j].y, currentColumn[j].z, 
                          faultIndex + 1);
                Point3D p2(nextColumn[j].x, nextColumn[j].y, nextColumn[j].z, 
                          faultIndex + 1);
                Point3D p3(nextColumn[j+1].x, nextColumn[j+1].y, nextColumn[j+1].z, 
                          faultIndex + 1);
                Point3D p4(currentColumn[j+1].x, currentColumn[j+1].y, currentColumn[j+1].z, 
                          faultIndex + 1);
                
                addFaultQuad(p1, p2, p3, p4);
                quadCount++;
                
                logFile << "创建四边形 " << quadCount << ": (" 
                       << p1.x << "," << p1.y << "," << p1.z << ") - "
                       << "(" << p2.x << "," << p2.y << "," << p2.z << ") - "
                       << "(" << p3.x << "," << p3.y << "," << p3.z << ") - "
                       << "(" << p4.x << "," << p4.y << "," << p4.z << ")" << std::endl;
            }
        }
        
        logFile << "断层 " << fault.name << " 共创建 " << quadCount << " 个四边形" << std::endl;
    }

    logFile << "\n断层面片统计：" << std::endl;
    logFile << "总四边形数: " << faultQuads.size() << std::endl;
    logFile << "总顶点数: " << vertices.size() << std::endl;

    logFile.close();
}

void SurfaceExtractor::extractSurfaces() {
    // 首先创建断层面片
    createFaultSurfaces();
    
    // 然后处理地层界面
    const auto& gridSize = gridModel.getGridSize();
    const auto& coordData = gridModel.getCoordData();
    const auto& zcornData = gridModel.getZcornData();
    
    std::vector<std::vector<std::array<double, 6>>> coord_reshape;
    coord_reshape.resize(gridSize.Ny + 1);
    for(int j = 0; j <= gridSize.Ny; j++) {
        coord_reshape[j].resize(gridSize.Nx + 1);
        for(int i = 0; i <= gridSize.Nx; i++) {
            int idx = (i + j * (gridSize.Nx + 1)) * 6;
            for(int k = 0; k < 6; k++) {
                coord_reshape[j][i][k] = coordData[idx + k];
            }
        }
    }

    for (int k = 0; k <= gridSize.Nz; k++) {
        size_t baseIndex = k * 8 * gridSize.Nx * gridSize.Ny;
        if (k > 0) {
            baseIndex -= 4 * gridSize.Nx * gridSize.Ny;
        }

        int nodeType = (k < 10) ? 11 + k : 0;

        for (int j = 0; j < gridSize.Ny; j++) {
            for (int i = 0; i < gridSize.Nx; i++) {
                size_t idx = baseIndex + 2*i + 2*j*gridSize.Nx*2;
                
                double z_values[4] = {
                    zcornData[idx],
                    zcornData[idx + 1],
                    zcornData[idx + 2*gridSize.Nx],
                    zcornData[idx + 2*gridSize.Nx + 1]
                };

                const auto& coord00 = coord_reshape[j][i];
                const auto& coord01 = coord_reshape[j][i+1];
                const auto& coord10 = coord_reshape[j+1][i];
                const auto& coord11 = coord_reshape[j+1][i+1];

                Point3D vertices[4];
                
                for(int v = 0; v < 4; v++) {
                    const auto& coord = (v == 0) ? coord00 : (v == 1) ? coord01 : (v == 2) ? coord10 : coord11;
                    double t = (z_values[v] - coord[2]) / (coord[5] - coord[2]);
                    t = std::max(0.0, std::min(1.0, t));
                    double x = coord[0] + t * (coord[3] - coord[0]);
                    double y = coord[1] + t * (coord[4] - coord[1]);
                    vertices[v] = Point3D(x, y, z_values[v], nodeType);
                }

                addStrataQuad(vertices[0], vertices[1], vertices[3], vertices[2]);
            }
        }
    }
}

void SurfaceExtractor::addStrataQuad(const Point3D& p1, const Point3D& p2, const Point3D& p3, const Point3D& p4) {
    size_t v1 = addUniqueVertex(p1);
    size_t v2 = addUniqueVertex(p2);
    size_t v3 = addUniqueVertex(p3);
    size_t v4 = addUniqueVertex(p4);
    strataQuads.push_back({v1, v2, v3, v4});
}

void SurfaceExtractor::addFaultQuad(const Point3D& p1, const Point3D& p2, const Point3D& p3, const Point3D& p4) {
    size_t v1 = addUniqueVertex(p1);
    size_t v2 = addUniqueVertex(p2);
    size_t v3 = addUniqueVertex(p3);
    size_t v4 = addUniqueVertex(p4);
    faultQuads.push_back({v1, v2, v3, v4});
}

size_t SurfaceExtractor::addUniqueVertex(const Point3D& point) {
    auto it = std::find_if(vertices.begin(), vertices.end(),
        [&point](const Point3D& p) {
            const double epsilon = 1e-6;
            return std::abs(p.x - point.x) < epsilon &&
                   std::abs(p.y - point.y) < epsilon &&
                   std::abs(p.z - point.z) < epsilon;
        });

    if (it != vertices.end()) {
        if (point.type >= 1 && point.type <= 10) {
            it->type = point.type;
        }
        return std::distance(vertices.begin(), it);
    }

    vertices.push_back(point);
    return vertices.size() - 1;
}

